Coulomb matrix elements for the impact ionization process in nanocrystals: the 

envelope function approach 
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We propose a method for calculating Coulomb matrix elements between exciton and biexciton 
states in semiconductor nanocrystals based on the envelope function formalism. We show that 
such a calculation requires proper treatment of the Bloch parts of the carrier wave functions which, 
in the leading order, leads to spin selection rules identical to those holding for optical interband 
transitions. Compared to the usual (intraband) Coulomb couplings, the resulting matrix elements 
are additionally scaled by the ratio of the lattice constant to the nanocrystal radius. As a result, 
the Coulomb coupling between exciton and biexciton states scale as l/R 2 . We present also some 
statistical estimates of the distribution of the coupling magnitudes and energies of the coupled 
states The number of biexciton states coupled to exciton states form a certain energy range shows a 
power-law scaling with the ratio of the coupling magnitude to the energy separation. We estimate 
also the degree of mixing between exciton and biexciton states. The amount of biexciton admixture 
to exciton states at least 1 eV above the multiple exciton generation threshold can reach 80% but 
varies strongly with the nanocrystal size. 
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I. INTRODUCTION 

Limitations of the efficiency of the existing solar cells 
motivate continuous search for new technological solu- 
tions that might lead to more efficient photovoltaic con- 
version. One of the fundamental limitations on the ef- 
ficiency of the existing solar cells results from the fact 
that photons with energies higher than the energy gap of 
the semiconductor excite high energetic electrons in its 
conduction band. The excess energy of these charge car- 
riers is dissipated in a phonon relaxation processes and 
is therefore lost for photovoltaic conversion. 

Among systems that are investigated for possible use 
in solar energy conversion devices in order to overcome 
these losses are semiconductor nanocrystals (NCs). The 
increased efficiency of photovoltaic conversion in NC- 
based devices might result from multiple exciton gen- 
eration (MEG) due to impact ionization (inverse Auger) 
processes^. Such an effect consists in generation of two or 
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FIG. 1: Schematic representation of the impact ionization 
process in a nanocrystal. The state labels refer to labeling of 
the matrix elements derived in Sec. [Ill 



more electron-hole pairs by a single high energy photon 
and thus converts the excess above-bandgap energy into 
useful current. This process is enabled by Coulomb cou- 
pling between single-pair (exciton, X) states and two-pair 
(biexciton, BX) states in a NC (or, in general, between 
states with n and n + 1 pairs), as shown in Fig. [T] It can 
result from the system dynamics after optical creation of 
a single electron- hole pair, during which an electron re- 
laxes within the conduction band transferring its energy 
to an intraband excitation (an inverse Auger process, in 
which the two situations in Fig. Q] are understood as "ini- 
tial" and "final"). However, very short time scales of the 
biexciton generation^ suggest that the process may be 
instantaneous and result from the mixing between the X 
and BX states. In this case, the original X states contain 
on the average more than one electron-hole pair, while 
the initially dark BX states become optically active. The 
MEG process can be thus understood as an excitation of 
a BX state mediated by a "virtual" X state to which it 
is coupled^ (the two states in Fig. [T] are then interpreted 
as the "virtual" state that is directly coupled to light 
and the final BX state that is optically excited via the 
Coulomb- induced mixing). Clearly, in this picture, the 
degree of mixing between the states and their spectral 
distributions (spectral densities) are crucial for the effi- 
ciency of the MEG process. Another possible process^ 
involves Coulomb coupling between empty NC and BX 
states and relies on an intermediate BX state. 

The MEG effect has been observed using a variety 
of experimental techniques in NCs formed of various 
narrow-gap semiconductors and under various excitation 
conditions^ - — , as well as in systems of coupled NCs^. 
In some experiments, quantum efficiencies as high as 
700% were observed^. On the other hand, some later 
experimental investigation showed much lower efficien- 
cies of the MEG process or even no traces of MEG at 
alli2r— . The subsequent discussion of the experimental 
factors involved in extracting the MEG efficiency from 
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experimental data^— showed that performing reliable 
experiments and correctly interpreting their results is not 
quite trivial. Uncontrolled effects present in the exper- 
iments, like photocharging of the NC core and charge 
trapping at the surface ligands, can indeed lead to a con- 
siderably overestimated result. In any case, however, 
enhanced photocurrent due to MEG has been observed 
in real NC-based solar cell o 20 ' 21 providing a direct proof 
of the usefulness of this process in solar energy conver- 
sion. Still, the experimental results and the controversies 
they arouse motivate theoretical work aimed at full un- 
derstanding of the MEG process in NCs. 

Theoretical description of the X and BX spectrum 
and the X-BX couplings is sought within various ap- 
proaches available for modeling semiconductor nanos- 
tructures. Atomistic approaches to this problem include 
density functional theor y 22 ' 23 , pseudopotential 2 ^— and 
tight binding22r— methods. However, the high compu- 
tational complexity of atomistic methods limits the size 
of tractable systems and forces one to radically truncate 
the basis of single particle functions when computing the 
properties of few-particle complexes. Therefore, one of- 
ten resorts to an envelope function approach based on 
the k-p theory 3 ^ - — , which was very useful in the past for 
understanding the fundamental electronic properties of 
NCs^i. Apart from establishing the spectral distribution 
of the X and BX levels, the central point of any theoreti- 
cal modeling of the MEG process is the evaluation of the 
X-BX Coulomb matrix element which can then be used, 
e.g, to find the spectral properties of correlated X-BX 
states 3 ^ or to study the system kinetics 3 ^" 3 ^^— . 

In the present paper, we present the calculations of 
the Coulomb matrix elements between X and BX states 
within an envelope function approach. This approach is 
much more explicit than the atomistic computations and 
yields relatively transparent formulas for the matrix el- 
ements that offer much more insight into the properties 
of the X-BX couplings. This allows us to point out that 
the correct treatment of the Bloch parts of the carrier 
wave functions (which was not made explicit in previous 
works) leads to selection rules for the newly created pair 
identical to those holding for optical transitions. In ad- 
dition, the interband character of the coupling reduces 
the matrix element (as compared to the usual electron- 
electron interaction) by a factor a/R, where a is on the 
order of the lattice constant and R is the radius of the 
NC, and leads to the 1/i? 2 dependence on the NC size. 
On the practical side, the envelope function approach 
leads to formulas that, for a spherical NC, can be eval- 
uated at very modest computational cost. This allows 
us to find coupled pairs of X-BX states in a very broad 
energy window and to study the distribution of the mag- 
nitudes of the matrix elements vs. the energies of the 
coupled configurations. Here, we restrict ourselves to the 
simplest example of an application of the proposed for- 
malism and calculate the matrix elements using a very 
simple model of wave functions in order to extract very 
general statistics on the X-BX couplings. The important 



conclusion from this study is that the X-BX state mixing 
is on the average dominated by relatively strongly cou- 
pled pairs of energetically close states. This means, on 
one hand, that admixture of the BX states to the X states 
(or conversely) cannot be treated perturbatively. On the 
other hand, it shows that truncating the set of states 
to an energy window around a given X (or BX) state is 
a well-defined and convergent approximation. Quantita- 
tively, the matrix elements found for an InAs NC are in 
the meV range, which is clearly very small compared to 
typical separations of single particle states in these struc- 
tures. Still, admixture of BX states to a given X state 
can reach several to a few tens per cent due to relatively 
large density of coupled BX states. 

The paper is organized as follows. In Sec. HH we de- 
rive general formulas for the Coulomb matrix element be- 
tween X and BX configurations in a spherical NC. Next, 
in Sec. IIIII we present a simple example of an applica- 
tion of this formalism. Sec. IIVI contains results on the 
statistical distribution of the X-BX couplings and an es- 
timate of the degree of X-BX mixing. Sec. [V] contains 
concluding discussion of the results and the outlook for 
their extension. 



II. MODEL AND GENERAL DESCRIPTION 

In this section we present a general derivation of the 
matrix element underlying the MEG process in a NC. 
The calculation will be performed in the envelope func- 
tion approximation. To simplify the notation we assume 
identical envelopes for all bands but generalization to ar- 
bitrary envelope wave functions is straightforward. 

A. General considerations 

We calculate the matrix elements between single- 
subband wave functions. The result can be generalized 
to derive the Coulomb couplings between states in more 
accurate models including band mixing. We consider a 
spherical NC of radius R with the dielectric constant e s , 
emerged in a host material with the dielectric constant 
£h- The electron-trion coupling responsible for the im- 
pact ionization process results from the Coulomb cou- 
pling between the states in conduction and valence bands, 
corresponding to an intraband transition of an electron 
accompanied by the creation of an additional electron- 
hole pair (inverse Auger process) as shown schematically 
in FiglTJ The relevant part of the Coulomb Hamiltonian 
is 

H - J,(ffiffaACT4) ffl t ffl t n . „ _i_H c 

1/1 1/2 1/3 1/4 AcTlCT2CT4 

where a aiUi , a\. v . are the annihilation and creation oper- 
ators for an electron in the conduction subband <r, and 
in the envelope state while a\ Ui , at are the annihila- 
tion and creation operators for an electron in the valence 
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subband A and in the envelope state v. The index z/, rep- 
resents the whole set of relevant quantum numbers. The 
graphical interpretation of this term is given in Fig. [1] 
We write the wave functions in the single-band envelope 
approximation in the form , J' (T1/ (r,s) = ip v (r)ua-(r, s), 
where ipv(r) is the envelope and Ua-(r,s) is the lattice- 
periodic Bloch part (s denotes the value of the electron 
spin). The position vector is written as r = r.; + £, 
where j"j denotes the center of the ith unit cell (u.c.) 
and £ lies in the first u.c. The integration over the 
whole NC volume is then expressed as a summation over 
all the unit cells and integration over one cell. We as- 
sume that the envelope functions change slowly, so that 
Vvfa + £) ~ Vvfa) and use the periodicity of the Bloch 
functions u a (ri+£) = ?v(£). The matrix element is then 
given by 



„ „/ • ■ J u.c. J u.c. 



(1) 



x< faX s')C 2 (n)< 2 (£, 8)U(n + 1 rj + £') 
x^fa)«A(£, s)Vv 4 faX(£',s'). 

The two-particle interaction energy in a spherical NC is 
composed of the direct Coulomb interaction and the cou- 
pling via surface polarization due to dielectric disconti- 
nuity between the NC and the environment, U (r, r 1 ) — 
^dircctfa r') + U po i(r, r'), where 



^direct far') = 

f/ P oi(r,r') = 



1 



4ire e s \r - r'\ 

„2 („ m l\k 



(rr') 

ak MfcTT Pfe ( C0S ^)- 



47re e s R 2k+ 



Here cos 7 = r ■ r' / (rr 1 ), are Legendre polynomials, 
Xk = (k + l)(e — l)/(ke + k + 1) with e = e s /eh, and we 
take into account only the two-particle part of the surface 
polarization term. 

In consequence, the Coulomb matrix element splits 
into the corresponding two contributions 



I,(CTlcr2 Atr 4 ) _ t (cricr2Aer4,dir) , l(<ti<T2Act4,Po1) 



B. Direct Coulomb coupling 



(2) 



For the direct Coulomb term, which is singular at 
r = r', we split the summation into % = j (short-range 
contribution) and i ^ j (long range contribution). For 
the former, we use the expansion^ 



1 00 n 1 
t =Vi<. v 



1=0 ^> m=-l 



47T 

21 + 1 



Y^(9'<l>')Y lm (0,<l>), (3) 



where £< = min(£, £'), £> = max((, £'), 8, <f> and 9, <p' are 
the spherical coordinates of the vectors £ and £', respec- 
tively, and Yi rn (9, cj>) are spherical harmonics. Since both 
<7i and 04 correspond to s-type conduction band states 
only the term with I = is non-zero in the integral over 



£' when Eq. (J3J) is substituted to Eq. (JTJ. This term, 
however, yields a vanishing integral over £, as A and 03 
correspond to states with p and s symmetries, respec- 
tively. Hence, the short-range contribution vanishes. 

In the long range term, we expand |r< + £ — rj — 
to the linear order in (£ — £')• The summation over i,j 
in Eq. ([T]) is then replaced by integration, where one for- 
mally excludes a small volume around r = r' (which is 
represented by a prime over an integral)^. Thus one 
finds in the leading order 



p 2 1 

. O lC r 2 Acr4,dir) _ e J_ 

»iw>* 47re e s ^ V 2 



(4) 



xj d 3 rj d 3 rV; i (r')C 2 (^)V^— ^ 3 fa^V 4 fa) 

x / d 3 r<(£',s'K 4 (M / d 3 e< 2 (^,s)|u A (|, S ), 

where V is the volume of the unit cell. The zeroth order 
term as well as the term containing £' vanish because the 
functions w^j and Ug- 4 both have s-type symmetry. The 
integration over r' can be extended onto the whole space 
because the singularity of V(l/r) is integrable in three 
dimensions. Using the orthogonality of Bloch functions 
and the definition 



E 



(5) 



one then finds 



L(<Ti<T 2 AtT4,dir) 

" - i/ii/ 2 i> 3 i/ 4 



« e' / ,3 g ' a ^A 

d q — - — (6) 



(2tt)3 e e s 7 ~* q 2 

X J~ i> 1 v i \q)J~ : 1/31/2(9)^0-1 Oil 

where we used the identity 

1 _*_ f d 3 q ^l (r'-r)^ 



\r-r'\ 2tt 2 
and the form-factors are defined as 



T vv ,{q) = J d 3 rr v (r)e iq - r ^'(r). 



(7) 



Note that, according to Eq. ([6]), the leading order term 
in the matrix element responsible for the MEG (impact 
ionization) process involves V(r — r')^ 1 (as opposed to 
just (r — r') _1 in an intraband matrix element), hence 
the resulting quantity is proportional to 1/i? 2 . This is 
a consequence of the orthogonality of the Bloch func- 
tions that leads to the appearance of the bulk interband 
matrix element of the position operator, a„ 2 \, (which 
obviously has the dimension of length) in the interband 
Coulomb term. Note that the resulting selection rules 
for the impact ionization process are the same as for the 
dipole-allowed optical transitions. 

In a crystal with zinc-blende structure, the topmost 
valence band corresponds to the 3/2 band angular mo- 
mentum. The non-zero matrix elements a a \ between 
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the four valence bands with A = ±1/2, ±3/2 and the two 
conduction bands with a = ±1/2 are 



a_i_j_ _i_3 = \/3ct^i +i 

^2^2 "T 2 1 



a+i +i 




(8b) 



where ao is a bulk material constant. Hence, we find 



q ■ a<j 2 .\ = acr 2 \qa d -^-Yi,Am(#, tp), (9) 
where Am = A — a 2 and the non-zero coefficients are 



_i_3 = 1, a^-i 

a= 2' :c 2 +2' =E 2 



n/3' 



where 



L(cricr 2 Acr4,pol-j) 
^1^2^3^4 



(10) 

^'^^rE^ / dVl / ^^(^(n) 

x^^ 1 E 

J ;</e m=—l 

l—k-\-j odd 
x !/>„ (6» 2 , 02 ) VV 3 (»*l)^4 ( r 2 )5o-! , CT4 

for j = 1,2. 

As previously, using Eq. (|8a| and Eq. (|8b[). one finds 



(11) 



From Eq. (dJ), again the general scaling of the matrix 
element as 1/R 2 follows. 



C. Surface polarization contribution 

The surface polarization term is smooth, therefore sep- 
aration into short-range and long-range parts is not nec- 
essary. Again, we expand the potential into the Taylor 
series in £, 

C/poi(n +€,rj + £') » U po i(ri,rj) +£-'ViU po i(r i ,r j ) 
+ £'-VjU po i{r i ,r j ) + ..., 

and keep only the lowest order non- vanishing term. The 
first and third terms in the above expression lead to van- 
ishing integrals in Eq. Jl]) because of the orthogonality of 
the Bloch functions for different bands. Using the prop- 
erties of Legendre polynomials and the addition theorem 
for spherical harmonics^ we find for the second term 



€ ■ ViUpoiir^Tj) = 



e e s R^ R 2k 



n 



E E YhnVMYlmWi'h) 
l<k m— — l 



E E y . 



7 

l — k odd 



where (n, 6i,<f>i) are the spherical coordinates of the vec- 
tor 7*j. One substitutes the resulting expression into 
Eq. ((T|), performs the integration over the Bloch func- 
tions according to Eq. ([5|) and changes the summation 
over i,j into integration as previously. As a result, the 
surface polarization contribution to the matrix element 
is 

^(cio^Ao^jpol) _ L (o^O^Ao^, pol—1) I L (CTiCT 2 Acr4,pol— 2) 



D. Coulomb coupling between two- and 
four-particle configurations 

Assuming that the NC is spherical and neglecting the 
spin-orbit coupling, the total spins of electrons and holes 
are separately good quantum numbers. The bright states 
involving the hole in a state with ±3/2 angular momen- 
tum will be denoted by \a/3 an d \ a P 4-lt)> where a 
and /3 denote the electron and hole states (representing 
the relevant sets of quantum numbers) and the arrows 
represent the values of the projection of the band angu- 
lar momentum of the electron and hole on the selected 
quantization axis. The states with dark spin configura- 
tions as well as the states involving a hole with ±1/2 band 
angular momentum (represented by f, J.) are denoted in 
an analogous way. Note that, e.g., ft denotes hole spin 
±3/2, which results from a transition from the A = —3/2 
valence band. 

The four-particle states are labeled by the electron spin 
configuration (S,T±^) for singlet and the three triplet 
states, respectively and by the hole spin configuration. 
If only the topmost valence band (j = 3/2) is included 
than the two-hole configurations can be classified in the 
following way, which is convenient for our purpose: The 
states with both holes in ±3/2 or both holes in ±1/2 
states are combined in singlet-like and triplet-like con- 
figurations S 3 / 2 ,T^ and S" 1 / 2 ,^. For instance, in 
terms of the hole creation operators h A , the two- hole 
state with S 3 ^ 2 and T^ 2 spin configurations is |/i/i'Sh) = 

(w/^X^frfyU ± fyU^U)l°>> where M > A*' (equal- 
ity allowed only in the singlet configuration) and 77^/ = 
l/\/2 for fj, — /j', and = 1 otherwise. While these 
spin configurations are not total spin eigenstates, they 
have a definite parity under particle permutation, imply- 
ing also a definite (opposite) parity of the orbital wave 
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functions. Hence, the resulting states automatically di- 
agonalize the hole exchange interaction. The spin config- 
urations with one hole in a ±3/2 state and one hole in the 
±1/2 state are obtained simply by symmetrizing (S) or 
antisymmctrizing (A) the corresponding two-hole states 
with respect to the orbital wave functions. For instance, 
!////,%) = (w/V2)(^ t ^ 4 + ^ t ^ 4 )|0), // > H', 
and W,A n ) = (w/V2)(^ f ^, 4 - ^ t ^)|0), // > 
p! . The four-particle (biexciton) states are then labeled 
by |i/i/£ e ; ////'Eh), where v,v' denote electron states, 
//,//' are hole states, and E e ,Eh represent spin config- 
urations. 

The couplings between the two-particle and four- 
particle configurations, which are responsible for the im- 
pact ionization process, are then expressed in terms of 
the Coulomb matrix elements given by Eqs. © and 
(|10[) . For instance, the non-zero couplings between the 
four-particle states with both holes in the spin-±3/2 
states and the two-particle state \af3 fJJ.) are (denoting 
hx 6 z h = {w'Z e ,w'E h \H\apil}) 



ha a! 



SS 3 / 2 



(hi + h 2 ) + <V/3 ( h 3 + hi)] 



3/2 



ST. 



h ToS z/2 



TnT, 



0*0 



3/2 



-y- [5,,,p {-hi - h 2 ) + 5 pip (h 3 + hi)] 
%^ [8„p (hi - h 2 ) + <V/3 (h 3 - hi)] , 



[<W (~hi + h 2 ) + <V/3 (^3 - hi)] 



h 3/2 = Srf (— h'i + h' 2 ) + 8^0 (/ig - /I4) , 



h (tm) h 



where hi 

} h! 



h (tm) h 



"■vv'tia > hi 



,pt) 



K>'lZa ■ Since one 01 the holes is a spectator in the impact 
ionization process and its state is conserved there is no 
coupling between this two-particle state and any state 
with both holes in a ±1/2 spin state. The results for 
the state \aj3 J/fr) are obtained by flipping all the spins 
and the couplings for a two-particle state with a hole in 
a ±1/2 state are easily derived by exchanging the role of 
the ±3/2 and ±1/2 hole spins. 

With similar notation as above and A =j, I, the non- 
vanishing couplings between four-particle states with one 
hole in a ±1/2 state and one in a ±3/2 state and the same 
two-particle state are 



hsSxu = 



h SA* 



[<W (h'i + h'i) + 5^ (K + h'i)] 



hT Sx 



hT A>, 



hr+Sx 



l T+A^ 



2 



1 



[Stf (h'( + h'i) + 8^0 (h'i + h'i)] 

[<W (K-K) + 5,, p (h'i -h'i)]. 



^[8^(-K + K) + 5,, p (h'i -h'i) 



V2 
1 

71 



[<W (K 1 -K') + 8^,0 (h'i' -h'i')] 
[8^0 (-h'i' + h'i') + V (h'i' -h'i')] 



where h'( and h"' are defined as hi and h\, respectively, 
but with the third spin (upper) index replaced by A. 
Again, the results for the two-particle state | a/3 Tfl") are 
obtained by flipping all the spins. 

Thus, we have characterized the Coulomb couplings 
between two-particle and four-particle configurations 
within the envelope function approach. The matrix el- 
ements for transitions between various subbands can be 
used in a calculation of couplings in an arbitrary envelope 
function model, including the common one based on an 
8-band k ■ p Hamiltonian (obviously, the classification of 
the four-particle configuration must then be extended to 
account or the spin-orbit coupling). This rather complex 
task is beyond the scope of this paper. In the following 
section, we limit ourselves to the simplest application of 
the results obtained above to a single-band model. 



III. MATRIX ELEMENT FOR A SIMPLE 
MODEL OF NANOCRYSTAL WAVE FUNCTIONS 

In this section, as the simplest example of an applica- 
tion of the general envelope function formalism of Sec.llll 
we find the Coulomb matrix elements between the 2- 
and 4-particle states assuming simple single-band car- 
rier wave functions. While perhaps not quantitatively 
accurate, this calculation yields a useful estimate of the 
overall magnitude of the couplings and of the resulting 
degree of mixing between two- and four-particle configu- 
rations as well as some statistics of the coupling strengths 
that, taken globally, may be close to the actual ones. Be- 
low, we define the simple model of wave functions which 
is then used to implement the general findings of the pre- 
vious section. 



A. Model of the wave functions 

An InAs NC is modeled as a spherical potential well 
of radius R with infinite potential walls. The envelope 
wave function is given by 



i>u(r) 



1 



R 3 /- 



■N nl Yi m (6,cl>)ji(xi n r/R), 



(12) 



where we write explicitly v = (nlm), xi n is the nth 
zero of the spherical Bessel function ji and N n i — 
\2/\ji + x(xi n )]. In our model, we neglect band mixing 
and include the Coulomb interaction between electrons 
and holes in the lowest order only. 

The corresponding energy levels for electrons and holes 



(e,h) 
nl 



are E, 

effective mass m? 



h 2 x 2 n /(2m* cYl R 2 ). We take the heavy hole 



0.35mo, where too is the free 
electron mass, and use the implicit formula for the 
energy-dependent electron mass in the decoupled bands 
approximation^ m* = mo [a + Ep/(E g + E n i)] , where 
E P = 22.2 eV, E g = 0.418 eV is the bulk band gap, 
and the parameter a — 0.77 accounts for the coupling to 
higher bands. Using the implicit formula for the electron 
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and a constant effective mass for the hole is motivated 
by the small effective mass of the former which leads to 
larger kinetic energies as compared to the hole. 

In the energy levels of the few-particle states, we in- 
clude the lowest order corrections due to Coulomb inter- 
actions, including the surface polarization terms^S. In 
the energies of the four-particle states also exchange in- 
teractions are taken into account. 



B. Direct Coulomb coupling 

In Eq. ([7]), we substitute the wave functions from 
Eq. (fT2]) and use the expansion 40 



I in 

where (r, 8, (f>) and (q, i9, tp) are the spherical coordinates 
of the vectors r and q respectively. As a result we get 

i+i' 

FvAq) = M-i) m - m ' f 'fwv (qR) 

i"=\i-i'\ 



xG ™,™-™,rn Yiiim _ mi ^^ ^ 



where 



fu"v( u ) = NimNi'm' / dyy j v >{uy)ji{xi n y)j v (x Vn >y) 
Jo 



and 



Gram, l 
U'l" 



2-rr p 7r 



d<j) / d9 sin 8 
o Jo 

xY l * m (d,c/>)Y l ' m '(0,tt>)Y l "m"(9>< 



(14) 



are Gaunt coefficients. 

Using Eqs. ©, ©, and fH) we have 

u (CTia-2ACT4,dir) _ 

""V 1 7/ 2 V3 Ui 

4e 2 a 



3ire e B R z 

h+U h+h r x 

l-l'+l 



E E 

i=\h-i 4 \i>=\i 2 -i 3 \ 



pfni ,mi -m4,fii4 ^3 ,7713 — 7712 ,7772 /~irn\ —7774 ,±1.7773 — 7772 
XLr ii,M4 U h,l',h ' 

where Am = A — 02 ■ The 1 / R 2 dependence of the matrix 
element is explicit in this result. Note also that I — V must 
be odd for the Gaunt coefficients to be non-zero so that 
the matrix elements are real. 



C. Surface polarization contribution 

In the surface polarization- related term [Eq. (fTU|) ]. we 
use Eq. (1111) and substitute the wave functions from 



Eq. (TT2|) . To reduce the product of four harmonics, we 
expand one pair into the Gaunt series, 
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The resulting integrals are performed using Eq. (fT4"|) . As 
a result, one finds 
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and Am is defined as previously. Again, these results 
show an explicit 1/i? 2 dependence. 



IV. RESULTS 

In this section, we present results of calculations per- 
formed within the simple single-band model presented 
above. We focus on the general statistical distribution of 
the coupling strengths between optically active (bright) 
single-exciton (X) and biexciton (BX) states vs. the en- 
ergy differences between the two coupled states. This 
allows us to estimate the degree of mixing between the 
bright X and BX states. The calculations are performed 
for a single spherical NC with the radius R = Rq — 3 nm 
or for a slightly inhomogeneous ensemble of NCs with 
the sizes given by the Gaussian distribution of their radii 
with the mean Rq and the full width at half maximum 
(FWHM) of a = 0.3 nm. We take a = 1.9 nm, as esti- 
mated from the InAs parameters routinely used in k ■ p 
calculations^ 2 -. For definiteness, we focus on the bright X 
states with the hole with spin projection +3/2 that lie in 
the energy range below 5 eV. This set contains 53 states 
in a NC of 3 nm radius (the spectrum is showed in the 
inset in Fig. [2(c)). The lowest BX state, which sets the 
energetical onset of the MEG process, is at 2.6 eV. 



7 



10 

5 

-5 

-10 

8 
6 
4 
2 

-2 
-4 
-6 
-8 

10 

5 

-5 
-10 

6 
4 
2 

-2 
-4 
-6 
10 

5 



-5 

-10 



A = +3/2, 4m = ±1 ■ 
A = ±1/2, 4m = ±1 □ 
A = +1/2, 4m = x 



;xx 



(a) (2,1,-1) 



. ■ 



* 



(b) (2,2,1) 



E (eV) 23 
30 22 J5- 
14 




(c) (3,0,0) 



(d) (2,3,0) 




! "1 i:l!fPl 



(e) all states 



-400 -200 200 

AE (meV) 



400 



FIG. 2: (a-c) The matrix elements h between X and BX 
states vs. the energy distance between these states for se- 
lected X states, as shown in the labels, for a NC with the 
radius R — 3 nm. The point style encodes three groups of 
spin configurations of the BX state as shown in the panel (a) 
(see also the text), (d) The matrix elements and X-BX en- 
ergy distances for all the X sates below 5 eV. Inset in (c): the 
spectrum of the X states (the digits show the n and I values) . 



In order to characterize the typical X-BX Coulomb 
coupling strengths and the distribution of energies of 
the coupled BX configurations, in Fig. [2ja-c) we graphi- 
cally represent these couplings for three selected bright X 
states (indicated by the values of the quantum numbers 
(n, I, m), identical for the electron and the hole, shown in 
the label of each panel). Couplings to BX states within 
the energy interval of ±550 meV around the energy of a 
given X state are shown. Each symbol corresponds to a 
single BX state coupled to a given X state and its posi- 
tion shows the energy distance from the X state and the 
value of the Coulomb matrix element coupling the X and 
BX states. The BX states are divided into three groups 
according to their spin configurations: 5 configurations 
with the spin of the newly created hole A = ±3/2, 6 
configurations with A = ±1/2 and the envelope angular 
momentum change Am = ±1, and 6 configurations with 
A = ±1/2 and Am = ±1. These three groups are coded 
into the symbol styles, as shown in the key inserted in 
Fig.rJIa). 

One can see that typical values of the X-BX couplings 
are up to several meV but most of them are at most on 
the order of 1 meV (we have found a small number of 
stronger couplings, even over 30 meV, but only between 
energetically very distant states). There is a clear pat- 
tern in the spectral distribution of the coupled BX states 
which results from the shell structure of the NC spec- 
trum. An important feature is the growing number of 
the coupled BX states which is consistent with the rapid 
growth of the overall density of BX states with increasing 
energy. 

The features observed in the case of the three selected 
X states shown in Fig. [2j a-c) are confirmed by the analy- 
sis of the combined distribution for all the X states with 
energies below 5 eV, shown in Fig. [U[d). Altogether, 
there are almost 7200 BX states coupled to the 53 X 
states in this energy range, which is still only a tiny frac- 
tion of the total number of the BX states in this energy 
interval. 

The results discussed above suggest that the growing 
number of remote BX states that are sufficiently strongly 
coupled to a given X state can considerably contribute to 
the mixing between X and BX states. Hence, a picture 
based on BX states in vicinity of the X state (an "en- 
ergy window") may be misleading. In order to achieve 
a more complete picture we have found all the states in 
the energy range ±4 eV from each of the X states under 
consideration. This allows us to account for all the BX 
states for which the ratio of the X-BX matrix element h 
to the energy distance AE is greater than 0.01, except 
for the uncommon cases of \h\ > 40 meV (we have not 
found a single instance of such a large coupling). The 
ratio \h/AE\ is an important parameter as it determines 
(within the range of applicability of the perturbation the- 
ory) the admixture of the BX state into the X state which 
is crucial for the efficiency of the MEG process. In Fig. El 
we plot this ratio as a function of the energy of the X 
state involved. For more clarity, the results are divided 
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FIG. 3: The ratio of the X-BX matrix element h to the energy 
distance AE between the coupled states plotted as a function 
of the energy of the X state involved for a NC with the radius 
R = 3 nm. The label (n, I) indicates the quantum numbers 
and the label a/b shows the numbers of coupled BX states 
with \h/AE\ > 0.01 (a) and \h/AE\ > 0.01 (6). 



into three groups according to the spin configurations 
of the BX state, as previously. Each "stack" of points 
in this figure corresponds to a set of states with fixed 
values of the quantum numbers n and I, as denoted in 
the figure (with the exception of the states (3, 0, 0) and 
(1,5, m) which are accidentally almost degenerate for this 
NC size). As expected, the number of coupled BX states 
tends to grow with the energy of the X state which can 
again be attributed to the growing density of states of the 
BX states. As one can see, in spite of the enormous num- 
ber of BX states in the energy range taken into account 
the number of states with \h/AE\ > 0.01 can only reach 
several tens for a single X state (note that the numbers 
in the figure give the total number of BX states coupled 
to all the 21 + 1 X states). Moreover, only for a small 
fraction of these states one finds \h/AE\ > 0.1 (which 
corresponds to an admixture of BX state above 1%). For 
the specific NC size unser study only in one case a close 
resonance between an X state and a coupled BX state 
was found that resulted in the \h/AE\ ratio exceeding 1. 
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FIG. 4: Histograms of the number of coupled BX states as a 
function of the ratio q = \h/AE\. The color coding represents 
three groups of spin configurations of the BX states and the 
dotted lines show a 1/q 2 dependence, (a) A single NC with 
R = 3 nm. (b) An ensemble of NCs with the mean radius 
Ro = 3 nm and the FWHM of the radius distribution of 
0.3 nm. 



Comparing the small number of states with \h/AE\ > 
0.01 found in the ±0.55 eV energy window (as shown in 
Fig. [2]) with the much larger numbers found in the broad 
energy range (Fig. one might conclude that remote 
states play an important role in the X-BX state mix- 
ing and, in consequence, in the MEG efficiency. From 
the point of view of theoretical modeling, this opens the 
critical question whether the actual situation can be rea- 
sonably approximated by a model that takes into account 
only BX states in a certain, sufficiently large spectral win- 
dow around a given X state or, in other words, whether 
the approximate results converge sufficiently fast when 
extending the spectral window. In order to approach the 
answer to this question we have estimated the statistical 
distribution of the values of q = \h/AE\ based on our 
simplified NC model. The result for a single NC with 
R = 3 nm is shown in Fig. |4ja). In order to obtain a 
single characteristics we present a joint distribution of 
the \h/AE\ ratios for all the X states below 5 eV, di- 
vided again in the three groups of spin configurations. 
As could be expected, the number of weakly coupled BX 
states is the highest, while the number of states cou- 
pled by larger matrix elements decreases quickly with the 
growing strength of the coupling (note that the number 
of BX states in the lowest sector of q, not shown in the 
figure, is formally infinite). The same tendency is seen 
in Fig. 2]Jb) where an analogous distribution is shown in 
the same way for an ensemble of NC in which we have 
assumed a Gaussian distribution of the radii with the 
mean Rq = 3 nm and the FWHM equal to 0.3 nm. In 
both dependence of the form 1/q 2 is found, which 

is rigorously confirmed by the logarithmic plot in the in- 
set to Fig. |4|b) , where this power law dependence is seen 
to be maintained over a surprisingly broad range of the 
q values. 

The fact that the exponent of this power law distribu- 
tion is equal to 2 is remarkable: According to the per- 
turbation theory (which is valid for small values of q), 
the admixture of a single BX state to a given X state is 
equal to q 2 . Hence, if there are N q states with a certain 
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0.3 nm. 



value of q then their joint contribution is equal to N q q 2 . 
The scaling N q ~ 1/q 2 means that the growing number 
of weakly coupling states exactly compensates the de- 
creasing magnitude of the matrix element so that, on the 
average, BX states with all the coupling strengths con- 
tribute equally. This has the immediate consequence that 
discarding the part of the coupled states with q < g m ; n 
generates a computational error that decreases propor- 
tionally to <7mi n - This means that, in principle, a numer- 
ical computation can be performed with arbitrary accu- 
racy based on a properly selected set of BX states. In 
fact, since the values of the matrix elements seem to be 
bounded, this implies convergence in terms of the width 
of the energy window assumed. This property, that holds 
only for a power law exponent below 3, is very desir- 
able and usually implicitly taken for granted in numeri- 
cal computation o 24 i 25 ' 27 i 35 ' 36 i 43 but it seems by no means 
obvious. 

A direct confirmation of this "homogeneous contribu- 
tion" property resulting from the 1/q 2 scaling is pre- 
sented in Fig. [SJ Here each histogram bar is multiplied 
by q 2 . As a result, a remarkably constant distribution is 
obtained, which is particularly visible in the inset, where 
the cumulative distribution is shown, which has a linear 
form across a very wide range of values. While this flat 
distribution is an interesting feature, it should be noted 
that the value of q 2 N q has no direct physical meaning at 
high q where the perturbation theory is not applicable. 
In fact, this part of the distribution is due to very few 
BX states that come to resonance with some of the X 
states for a certain value of R. At resonance, the value 
of q is infinite but the maximum admixture in the case 
of just two resonant states is 1/2. 

Finally, let us estimate the degree of mixing between 
the X and BX states. In Fig. |6l we show the BX admix- 
ture to the selected X states considered here as a function 
of the NC radius. This is obtained separately for each X 
state by diagonalizing the Hamiltonian including only the 
X state in question and all the BX sates directly coupled 
to this X state by matrix elements for which q > 0.01. 
For each R, the eigenstate 'J'o with the highest X con- 
tribution is found and the BX admixture is determined 
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FIG. 6: Biexciton admixtures to selected exciton states. 



as 

j Pbx = EK s ^I*o)I 2 = i-K^|*o)| 2 , 

i 

where \X) is the X state and \BXi) are all the coupled 
BX states. We show the results for one selected state 
out of each subspace with fixed n and I, except for the 
states (2,2, m) which are presented in full to sec that the 
admixture to states that differ only by the value of m is 
very similar (which is due to the fact that these states 
are almost degenerate). 

As expected, the low energy states (1,0,0) and 
(1, l,m), lying below the theoretical MEG threshold of 
2.6 eV, have the BX admixtures well below 1%. However, 
the same is true for states lying up to 1 eV above the 
threshold. This property results from the selection rules 
that hold in the MEG process, as discussed in Sec.|TTJ An 
X state needs to be located in the vicinity of a BX state to 
which it is coupled, which becomes likely only when the 
density of BX states becomes sufficiently high. On the 
other hand, the states with energies above 4 eV are typ- 
ically much more strongly mixed with BX states. Even 
in these case, in spite of hundreds of BX states coupled 
to each such X state, the admixture only in some cases 
approaches 80% for a certain NC radius and for most 
(but not all) states drops down as R increases. Note 
that the admixture strongly depends on the NC size and 
in the case of the state (1, 5, 0) shows strong oscillations. 
In general, while the X states are clearly not completely 
dissolved in the densely distributed BX states, the degree 
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of mixing between X and BX states becomes considerable 
at energies a few eV above the MEG threshold which may 
lead to the MEG efficiency of a few tens per cent. 

V. DISCUSSION, CONCLUSIONS AND 
OUTLOOK 

The main result presented in this paper is the method 
for calculating Coulomb matrix elements between exciton 
and biexciton states in semiconductor nanocrystals based 
on the envelope function formalism. We have shown that 
such a calculation requires proper treatment of the Bloch 
parts of the carrier wave functions which, in the leading 
order, leads to spin selection rules identical to those hold- 
ing for optical interband transitions (however, the rules 
for envelope states are different). Moreover, the result- 
ing matrix elements are additionally scaled by the ratio 
of the lattice constant to the NC radius, as compared to 
the usual (intraband) Coulomb couplings. As a result, 
the Coulomb coupling between X and BX states scale as 
1/R 2 . 

Once the matrix elements between single-band states 
are known they can be used for calculating X-BX cou- 
plings using more exact carrier states found by diago- 
nalizing the 8-band Kane (fe • p) Hamiltonian. This ap- 
proach has been found to correctly reproduce the NC 
spectrum down to very small sizes^ and, when combined 
with the results presented here, can provide quantita- 
tively accurate description of the MEG process in NCs. 
For more accuracy, Coulomb correlations, in particular 
in BX states^ 5 ., could also be included. While atomistic 
models may offer more accurate single particle wave func- 
tions and allow one to include more system features (e.g., 
surface defects^ 5 -), their high computational costs limits 
the extent to which few-particle states and couplings be- 
tween them can be treated. Typically, when following 
an atomistic approach, one is forced to restrict the cal- 
culations to an energy window around a given X or BX 
state and to truncate the basis of BX states when simu- 
lating the system dynamics 35 . From this point of view, 
the method based on the envelope function formalism 
may offer a complementary approach to the trade-off be- 
tween the accuracy of single-particle states and the re- 
liability of few-particle modeling in which the accuracy 



of the atomistic models is sacrificed in favor of lowering 
the computational effort, which offers considerably more 
flexibility on the subsequent stages of theoretical anal- 
ysis, including the system dynamic s 35 ' 36 and dissipative 
evolutio n 31 ' 32 ' 37 ' 38 , where finding the X-BX coupling is 
the essential prerequisite for further modeling. An addi- 
tional benefit is the transparent nature of our envelope 
function method, which offers mostly analytical treat- 
ment and does not rely on large computational resources 
or dedicated software, hence can easily be employed by 
a broad community of researchers. 

While performing full multi-band k ■ p calculations is 
beyond the scope of this work, we have presented some 
preliminary estimates of the statistical distribution of 
the coupling magnitudes and the energies of the coupled 
states using a very simple single-band envelope function 
approach. Such overall statistical conclusions are likely 
to be valid even if the underlying characteristics of indi- 
vidual states are not absolutely accurate. We have shown 
that the number of BX states coupled to X states form 
a certain energy range scales as inverse square of the ra- 
tio of the coupling magnitude to the energy separation. 
This scaling property is remarkable as it guarantees that 
the contribution of remote states is finite and control- 
lable, which justifies limiting the computation to an en- 
ergy window (no matter what computational method is 
chosen)^^^^^. 

Finally, we have estimated the degree of mixing be- 
tween X and BX states as a function of the NC size. Very 
small BX admixture to X states has been found below 
and within 1 eV above the MEG threshold. Much larger 
mixing, reaching 80%, appears for higher-energy states. 
The amount of BX admixture to this states varies quite 
strongly when the NC radius is changed by a fraction of 
a nanometer. This may suggest that modeling based on 
a single NC size may not be representative for average 
properties of an ensemble. 
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